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We present the outline of an IBM 7090 machine program for the Monte 
Carlo estimation of the percolation probability for a variety of space lattices. 
The under lying theory is briefly summarized . 

I. STATEMENT OF THE PROBLEM 

Percolation processes deal with the transmission of a "fluid" (dis- 
turbance, signal, etc.) through a ''medium" (material, region, etc.) 
against impediment by random irregularities situated in the medium.' 
This paper considers the case where the medium is a regular crystal 
lattice in two or three dimensions, consisting of "atoms" (the vertices 
or sites of the lattice) and "bonds" joining specified pairs of atoms. The 
next section will specify the structure of the lattice more completely. 
The fluid originates at one or more atoms of the lattice, called the source 
atoms, and flows from atom to atom along the connecting bonds. How- 
ever, each bond (independently of all other bonds) has a fixed proba- 
bility p of being able to transmit fluid and a probability q = 1 — p of 
being blocked: these randomly situated blocked bonds constitute the 
random impediments to the spread of the fluid. We write P N (p) for the 
probability that the fluid will reach (or "wet", as we shall say) more 
than N other atoms besides the source atoms; and the problem is to 
estimate P(p) = limbec P N (p). We do this by estimating P N (p) for a 
suitably large value of N: it turns out that N ~ 2000 is sufficient in 
many cases. The present paper describes the general organization of an 
IBM 7000 program for obtaining a .Monte Carlo estimate of P N (p). 
The numerical results appear elsewhere." 

t Oxford University, Oxford, England. 

DO!) 



910 THE BELL SYSTEM TECHNICAL JOURNAL, MAY 1962 

II. STRUCTURE OF THE LATTICE 

We consider below the three-dimensional problem in such a way that 
it contains as a special case the two-dimensional problem. A regular 
three-dimensional lattice consists of a number of fundamental cells, all 
identical apart from their position in space. Each cell is specified by 
integer coordinates (x, y, z) representing its position in space (x, y, z = 0, 
±1, ±2, • • • ). Each cell contains a finite number of atoms, limited in 
our IBM 7090 program to a maximum of eight atoms per cell, and 
denoted by A, B, • • ■ , //. 

In the problem, as originally posed in the previous section, some of 
the bonds may be one-way (i.e., only able to transmit fluid in a specified 
direction) while others may be two-way (i.e., able to transmit fluid in 
either direction). It can, however, be proved theoretically that P N (p) 
is unaltered if a two-way bond is replaced by two one-way bonds of 
opposite directions. It is therefore enough to consider the case where all 
bonds are one-way bonds, and we confine ourselves to this case hereafter. 

Let us write T(x, y, z) for the atom of type T(T = A, B, ■•■ , H) 
in cell (.r, y, z). The lattice is regular in the sense that, if there is a 
bond from T(x, y,z)\o T*(x*, y*, z*), then there is a bond from T(x + 
/, y + m, z + n) to T*(x* + /,//* + m, z* + n) for any I, m, n = 0, 
±1, ±2, • • • . Therefore we need only specify the terminal atoms reached 
by each bond from each atom of the cell (0, 0, 0). The program limits 
the number of bonds from a given atom to a maximum of 12. There is 
no restriction that distinct bonds from a given atom shall lead to distinct 
terminal atoms; and thus we may, if we wish, have two or more bonds 
in the same direction between a given pair of atoms. 

The machine receives information about the lattice structure from a 
series of input cards, having the format described below. To each atom 
of the cell (0, 0, 0) there is a "structure" card, followed perhaps by one 
or two "continuation" cards. The format of a structure card is: 

Columns || 5-13 I 15-23 I 24-26 I 27-35 I 36 I 37-45 I 46 | 47-55 1 56 I 57-65 I 70-72 

Contents l| STRUCTURE | T(+0+0+0) | T<?>- | T(xyz) \ , \ T{x U z) \ , \ T(xyz) \ , \ T(x„z) \ CTl 

In such a card, T (in columns 15, 27, 37, 47, 57) stands for one of the 
letters A, B, ■ • • , H (not necessarily the same in each case); and each 
of the symbols x, y, z is a signed integer. The entries in columns 5-13, 
24-20, 36, 46, 56, 70-72 are in BCD (binary coded decimal). Instead 
of the' whole word STRUCTURE in columns 5-13, the single letter S in 
column 5 will suffice. Columns 15-23 specify the atom of cell (0, 0, 0) 
from which bonds lead to the atoms appearing in columns 27-35, 37-45, 
47-55, and 57-65 respectively. Thus the structure and card provide for 
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Oil 



up to four bonds from the given atom of cell (0, 0, 0). If there are fewer 
than four bonds, some of the T(xyz) will be left blank. The next four 
(or fewer) bonds are similarly specified on a first continuation card 
with format: 



Columns 
Contents 



5-7 

CT\ 



27-65 
Same format as structure card 



70-72 
CT2 



The last four (or fewer 
with format: 



bonds appear on a second continuation card 



Columns 
Contents 



o-7 
CT2 



27-65 

Same format as structure card 



70 
S 



In both continuation cards, columns 5-7 and 70-72 are in BCD. If the 
atom specified in columns 15-23 of the structure card has four or fewer 
bonds from it, both continuation cards are omitted and we replace CT1 
in columns 70-72 by S in column 70 ( in BCD) ; if it has between 5 and 8 
(both inclusive) bonds from it, the second continuation card is omitted 
and we replace CT2 in columns 70-72 by S in column 70. The final 
structure on continuation card is followed by an "end" card with END 
(in BCD) in columns 5-7. 

For example, the simple cubic lattice with a pair of one-way bonds 
(one in each direction) between each pair of nearest-neighbor atoms is 
specified by: 

STRUCTURE --K+0+0+0) TV)- A(+l+0+0), .-1 (+0+1+0). A (+0+0+1), A (-1+0+0) CT\ 

CT\ -K+o+o+o) ro- a(+o-i+o*, .-K+o+o-i) s 

END 

Similarly, the tetrahedral lattice (diamond crystal) with a pair of one- 
way bonds in each direction between nearest neighbors is given by: 

STRUCTURE .41+0+0+0) T0- B(+0+0+0), fl(+l+0+0), £(+0+1+0), fi(+0+0+l) S 
S fi(+0+0+0) TO- .11+0+0+0). .41-1+0+0), .4 (+0-1+0), A (+0+0-1) s 

END 

The two-dimensional cases arise when 2 = identically on all cards. 



III. INPUT TO THE MACHINE 

The complete input to the machine consists of the program deck, 
followed by (i) an identification card, followed by (it) a parameter 
data card, followed by (Hi) a set of structure and continuation cards 
(as described above), followed by (iv) an end card (as described above). 

The format of the identification card is: 



Columns 
Contents 



5-15 
PERCOLATION 



18-25 
CRYSTAL 



26-66 

Xame of crystal 



68-72 

PUNCH 
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Those entries arc all in BCD. If PUNCH appears in columns 68-72, the 
output from the machine will appear on the printout and will also be 
punched onto cards; if PUNCH is omitted from columns 68-72, the 
output will appear on the printout but cards will not be punched. 
The format of the parameter data card is: 

Columns II 5-8 I 10-12 I 13 I 1-1-16 I 19-23 I 25-28 I 31-37 I 39-17 I 48 I 49-57 

Contents || RUNS | '" | - I *" I LINES \ "" \ SOURCES \ T(xyz) | , | T{xyz) 

Columns 5-8, 13, 19-23, 31-37, 48 are in BCD. Columns 10-12 contain 
a first run number; columns 14-16 contain a last run number: the 
purpose of these two run numbers will be described presently. All run 
numbers must be positive integers. Columns 25-28 contain the number 
N, appearing as a suffix in the desired function P N (p)- Tne maximum 
value of N permitted by the program is 7000. 

IV. GENERAL CONDUCT OF THE CALCULATION 

The calculation carried out by the machine is at first sight rather 
different from the problem posed in the first section of this paper. The 
change of formulation simplifies the calculation without affecting the 
ultimate numerical answer. 

The calculation consists of a number of separate runs. In each individ- 
ual run, independent random numbers rj are assigned to each bond of 
the lattice. Each 77 is uniformly distributed between and 1. This process 
of assigning 77-values to the bonds replaces the original process of block- 
ing bonds, so that the question of a bond being blocked or not does not 
arise in the reformulated process. Consider a connected path of bonds 
ft A 1 * " i A on tne lattice, where the (necessarily one-way) bond 
0i leads from the atom to which the preceding bond 0,_i led (i > 1 ). 
Let 771 , 772 , • • • , Vk be tin* 77-values of the respective bonds ft , /3 2 , • • • , 
j8 t . Define the 77-value of this path to be the minimum of 171 , m ,■■'»»»* ■ 
Next consider any atom other than a source atom. Define the 7- value of 
this atom to be the supremum of all path 77-values, taken over all paths 
which lead from some source atom to the given atom. Finally define c„ 
to be the largest number such that there are more than n atoms, other 
than source atoms, whose 7-values are equal to or greater than c„ . 

Now c„ for any fixed n is clearly a random variable, depending upon 
the several 77-values assigned to bonds of the lattice. P„(p), regarded as 
a function of p, is the cumulative distribution function of the random 
variable 1 — c„ . 

The machine is programmed to calculate c n for each run. Thus the 
set of all runs provides a sample of values of 1 - c„ , and the empiric 
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distribution of this sample can ho taken as an estimate of the required 
function P u (p). 

So far we have regarded n as a fixed integer. Actually, the output of 
the machine on any one run is a table of c„ as a function of n for all 
n ^ N, where N is the number set on the parameter data card. Thus the 
complete calculation provides estimates for P„(p) for all n = 1 , 2, • • • , 
N and all g p ^ 1 . 

Succosive runs are numbered serially Ri , R t -f 1, • ■ • , R., , where ft, 
is the first run number specified on the parameter data card, and R-, is 
the final run number specified on the parameter data card. Hence 
R-2 — R\ + 1 is the sample size for each empiric distribution. The run 
number R y + i of the (i ' + 1 )th run is used to trigger off the genera- 
tion of random numbers -q assigned to bonds in this run. Hence a run 
may be repeated for checking purposes by repeating the run number; 
but, if a fresh and independent sample of 7,, is desired, the parameter 
data card must specify a set of run numbers which does not overlap the 
set previously used. Since three decimal digits are available for run 
numbers, the maximum sample size is 999. A sample of size about 100 
is usually adequate. 

V. OUTPUT FROM THE MACHINE 

The printout from the machine, also punched onto cards if ordered 
on the identification card during input, is as follows. 

The printout begins with a copy of all input data (excluding the 
program deck). Thereafter follows a table, whose columns are headed: 

RUN N0. STATE c VALUE 

Each row of the table has an entry in each of these four columns. An 
entry under RUN is the current run number R being computed (R t ^ 
R ^ /?*)■ The entries under N0. and c VALUE are respectively n and 
c„ , tabulated for ^ n ^ N omitting any values of n such that c„ = 
r,,_i (0 ^ n ^ A/). Thus the only values of c„ printed are new values 
less than all preceding values in the run. (Clearly c„ is a non-increasing 
function n by virtue of its definition.) Such new values are indicated 
by the prefix NEW. The value of c A - is, on the other hand, prefixed by 
FINAL and also suffixed by an asterisk to help in reading the output 
quickly. The entry under STATE is either FINISHED or T(J) BE 
C0NT., according to whether the run is complete or not. The only other 
possible entry under STATE is INHIBITED: this is a safety device, 
to be explained later. It means that certain technical circumstances 
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Table I — Extract of Printout 



HUN 


N0. 


STATE 


c VALUE 


23 


518 


T0 BE C0NT. 


NEW c = .748 


23 


1000 


FINISHED 


FINAL c = .748* 


24 





T0 BE C0NT. 


NEW c = .906 


24 


1 


T0 BE C0NT. 


NEW c = .803 


24 


5 


70 BE C0NT. 


NEW c = .791 


24 


11 


T0 BE C0NT. 


NEW c = .783 


24 


17 


TO BE C0NT. 


NEW c = .770 


24 


19 


T0 BE C0NT. 


NEW o = .762 


24 


20 


70 BE C0NT. 


NEW c = .715 


24 


21 


70 BE C0NT. 


NEW c = .711 


24 


1000 


FINISHED 


FINAL c = .711* 



(whoso occurrence is extremely unlikely) have arisen to prevent com- 
pletion of the run. In an inhibited run all values of c„ in the printout 
are valid: all that has happened is that n has been prevented from rising 
beyond a certain value, at which instant the run is automatically dis- 
continued. 

The extract of the printout shown as Table I will help clarify matters. 
It gives the end of the twenty-third and the whole of the twenty-fourth 
run for a computation on the simple cubic crystal with N = 1000. To 
find the value of c„ for a value of n not printed in the Table, take the 
value of c for the largest n less than the required n. In the above example, 
d = 0.791 in run 24. 

The machine stores values of c as 9-bit binary decimals. Hence the 
rounding error in c is about 0.001 and not 0.0005. 

vi. outline of the program 

What has been said so far contains everything that a user of this 
program needs to know. What follows in the remainder of this paper 
is an explanation of how the machine carries out the program, and is 
intended for those who are interested in programming techniques. 

The two main entities in any given run are (i) a number denoted by 
c, and (ii) a so-called "wet list" of atoms. An atom is qualified for 
membership on the wet list if it is a source atom or if its 7-value is not 
less than c. Normally the value of c is held constant and the machine 
recruits new members of the wet list. However, if a stage is ever reached 
during a run where no further recruits can be found with the existing 
value of c, then the machine reduces c by an amount just sufficient to 
ensure the existence of at least one fresh recruit. The run begins with 
c = 1 and a wet list containing just the source atoms. Since a given 
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atom in a given run has a given y-value (depending upon the ^-values 
assigned to bonds in that run), reduction of c can never disqualify 
existing membership of the wet list. Successively recruited members of 
the wet list (other than source atoms) are numbered 0, 1, • ■ • . A little 
reflection will show that at the moment, when the member numbered n 
is added to the wet list, the current value of c must be c n . The run is 
terminated as soon as the wet list contains N + 1 members. Thus the 
machine output is simply the result of printing any freshly reduced 
value of c against the number of the next member to be added to the 
wet list. 

This procedure would be unworkable if the machine had to examine 
all atoms of the lattice for this qualification to belong to the wet list. 
What makes the procedure workable is the observation that there can 
exist, no qualified fresh recruits if the existing value of c exceeds the 
7j-values of all bonds, which lead from some atom of the wet list to some 
atom not in the wet list. Let us call such bonds the outgoing bonds of 
the wet list. If the wet list has an outgoing bond whose 77-vaIuc is at 
least c, then this bond leads to an atom which is qualified for membership 
of the wet list. If the value of c has to be reduced, the new value of c 
is the highest 17-value of all existing outgoing bonds. Hence, at any 
stage of the calculation the machine need only examine the 17- values 
of the outgoing bonds of the currently existing wet list. 

Information about the current status of the wet list resides in a block 
of registers in the core storage of the machine, with three registers 
(denoted by E, K + I, E + 2) allocated to each atom of the list. We 
write 

B = (x,y, z, T) = [11, 11,11,3J; 

/.; + 1 = (A') = [36]; (1) 

E + 2 = (./<,</>, t) = [9, 12, 15]; 

to indicate that E contains four different quantities (denoted by x, y, z, 
and T respectively) occupying 11, 11, 11, and 3 bits of the 3(5 bits 
available in a single register. The contents of E + 1 and E + 2 are 
exhibited in a similar manner. If we wish to emphasize that we are 
talking about the nth atom A n (n = 0, 1, ■ ■ , N) in the wet list (ex- 
cluding the source atoms of the wet list), we place n as a suffix to any 
of the above quantities: thus 

(E+ 2)„ = (^,,0„,7r„). (2) 
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The symbols have the following meanings: (x„ , y„ , z„) are the three 
integer coordinates of the cell containing the atom A„ , and T„ is the 
type of the atom A n (T n = A, B, ■ ■ • , H in the notation for crystal 
structure). X n is a 36-bit pseudo-random number, generated according 
to the multiplicative congruential recursive relation (low multiplication) 

X n s^X_ (mod2 38 ), (3) 

where is a permanently fixed odd integer, and the recursion is triggered 
from the initial state 

X_, = f(R) (4) 

where f(R) is a fixed function of R, the run number of the run under 
consideration. (R = Ri + j m the (j + 1 )th run). The meanings of 
\f/, <f>, and ir will be stated in a moment. 

The lattice structure allows us to have up to 12 bonds leading from 
each atom; and for the sake of exposition it is convenient to suppose 
in the first place that each atom has a full complement of 12 bonds 
from it. The atom A,, thus requires 12 77- values for its 12 bonds: we 
denote these by (i7„ a W 2 \ • ■ • , In™)- The quantity <f> consists of 12 
independent bits 0„ = (0„ (l, ,0,/' J) , • • • , 0n° 2) ) also corresponding to the 
12 bonds. If an atom has fewer than 12 bonds from it, we put the cor- 
responding <f> li> = to signalize the absence of a bond. For a bond 
actually present, we also put <f> U) = if it is known that this bond leads 
to some atom already in the wet list. In all other cases, ' =1. Thus 
all outgoing bonds from A„ have 0,/" - 1: the converse is not neces- 
sarily true, since we may have a bond, which exists and which leads to 
an atom already in the wet list (and is therefore not an outgoing bond), 
although at the current stage of the calculation we have not yet dis- 
covered that this bond leads to an atom in the wet list. Thus the quantity 
represents a state of current knowledge. We define 

1 / (1) , 0) d^l j C12)\ ( r\ 

ip„ = max (77,, </»„ , ■ • ■ , Vn 0» )• \ ,} i 

To validate this definition we require that each 77,/" be a 9-bit number. 
We achieve this by means of 

(X n ) - (vn ll) , Vn l2) ,vn 3 \ ?) = [9,9,9,9,] 

(»0-(^.*-.*. M .0-R«.9 > 9J (6) 

ifx n ) - u;^,%;v) - [9,9,9,9,] 

GMT.) - U, ,I,, W, , "U/ , V) = [9,9,9,9,] 
where ? denotes a 9-bit number which is not used (because the terminal 
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digits of pseudo-random numbers are too regularly distributed), and 
where the eongruential notation in ((}) indicates low multiplication as 
in (3). 

VII. PROGRESSIVE CONSTRUCTION OF THE WET LIST 

We are now able to describe recursively how the wet list is compiled. 
In what follows, it is important to remember that the wet list consists 
of both source atoms as well as recruited atoms in this list. 

.Suppose that the wet list is already partly compiled and that we have 
reached a stage at which the current c value has just been reduced to 
a new value. Starting at the beginning of the wet list, we successively 
scan each atom in the order of the list. For each atom scanned we ask 
first if its ^-value is less than c. If \p < c, we pass to the next atom on 
the list. If \p ^ r, we determine all values of i such that V'V" ^ c: 
these represent the only bonds which can lead to an atom at present 
qualified for membership of the wet list. Call such an atom a target 
atom of the scanned atom. Noting the cell coordinates and the type 
number of the scanned atom, we compute the cell coordinates and the 
type number of each target atom of the scanned atom. We then look 
through the wet list to see which of the target atoms do not belong to 
the wet list, and we add to the end of the wet list all target atoms not 
already on the wet list. With these new additions to the wet list, every 
bond from the scanned atom to a target atom leads to an atom of the 
wet list, and therefore we now set 1 " = for all values of i for which 
we had rf" <f> { " ^ c. Next we recompute \J/ from (5). Of course the new 
value of \J/ is less than r. Therefore in the scanning procedure we go to 
the next atom in the wet list. Ultimately the scanning procedure will 
reach the end of the wet list. At this stage, all \p in the wet list are less 
than c. We therefore reduce c to the largest ^ in the wet list, and restart 
the scanning from the beginning of the wet list. We continue this proce- 
dure until we have recruited N + 1 atoms to the wet list. 

To assist in computing the new value of c required at any reduction 
of c, we make a small modification of the foregoing procedure. We define 
a number c*, called the c-candidate. At the beginning of any scan, we 
set c* = 0. Before leaving any scanned atom and proceeding to the 
next one, we take max (\(/, c*) to be the new value of c*. Thus when we 
reach the end of the scan, c* equals the required new value of c. 

VIII. WET LIST SEARCH 

As described above, we have to search through the wet list to decide 
if a target atom is already in the list. To expedite this search we define 
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the modular typo of an atom A n in the wet list to be the least non- 
negative residue 

L = Xn + un + *» (mod 256), (7) 

and we dissect the wet list into 256 equivalence classes according to 
their modular type. A target atom can only be in the wet list if it is its 
own equivalence class, and hence it is .sufficient to search just one 
equivalence class in the wet list. 

Suppose that 4„ ( n , A n(2) , • • • , A n(k) are the atoms in the wet list 
currently belonging to the equivalence class to be searched. We define 
7r„ , appearing in (1), by 

* nlj) = n(j - 1). (8) 

Thus we can search the equivalence class backwards; for at any stage 
of the .search, the a-value for the atom currently examined gives us the 
address of the next atom to be examined. To start off this iteration, 
the core storage holds a table, called the //-table, with 256 entries pro- 
viding the values of n(k) for each equivalence class. When a new atom 
has to be added to the wet list, we must write up its ir-value: its v- value 
is simply the entry (for the appropriate equivalence class) that is in 
the //-table immediately before adding the new atom to the wet list; 
and, immediately after adding this new atom to the wet list into address 
n (/ ( . _|_ i ) f Ba y j W e enter n(k + 1) into the appropriate position of the 
//-table in place of n(k). 

IX. INHIBITION OF HUNS 

Each of the coordinates x, y, z of a lattice cell is represented by an 
11-bit integer, treated by the machine modulo 2048. Thus, effectively, 
the lattice lies on a four-dimensional torus instead of the required three- 
dimensional flat of four-dimensional Euclidean space. To remedy this 
defect, we cut the torus on each of the three two-dimensional flats defined 
by x = 0, y — 0, and z = respectively. We place the source atoms in, 
or in the immediate neighborhood of, the cell (1024, 1024, 1024); and 
we set an inhibition flag if any cell coordinate becomes zero modulo 2048. 

The run is allowed to proceed as before, after the inhibitor flag has 
been set, up to such time as the machine calls for a new value of c. At 
this instant, however, the inhibitor flag prevents the new value of c 
being set, and instead terminates the run with the printout comment c 
INHIBITED together with a point of the old value of c. 

The net effect of this procedure is to allow the fluid to pass across the 
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two-dimensional cuts in the torus, and even to complete circuits which 
are not nomotopic to zero. However, if such a circuit occurs, it may 
imply an unnecessary reduction in c, or an unnecessarily large reduction; 
and in that event the run must he terminated as a safety measure. 
Nevertheless it is very unlikely, with the values of N used, that inhibi- 
tion will be invoked; and in fact it has not been invoked on any run 
calculated to date. 

X. GENERATION OF PSEUDO-RANDOM NUMBERS 

There is no point in calculating pseudo-random numbers which are 
not going to be used. At the start of a run, all entries X, defined in ( 1 ), 
are set with a negative sign. When the machine comes to scan any atom 
in the wet list to look for possible outgoing bonds, it first asks if X is 
negative. If X is negative', it replaces A' by a positive pseudo-random 
number X„ , generated there and then by means of (3). Jf X is positive, 
the machine knows that a pseudo-random number has already been 
calculated for this atom, and it does not change A r n . Thus, as the prob- 
lem requires, each individual pseudo-random number remains fixed 
throughout a run. 

XI. DELTA AND PHI TABLES 

Consider the stage of the calculation when the machine is scanning the 
wet list and looking for the possible outgoing bonds from a particular 
atom A to the corresponding target atoms. Let (x, y, z, T) denote the 
cell coordinates and type number of .4 ; and suppose that we are con- 
sidering the ith bond from A (i = 1, 2, • • • , 12) as a possible outgoing 
bond to a target atom, whose cell coordinates and type number we 
denote by (x' t .//, z', 7"). Then 

/s" - K = (.,-', //, *', 7") - (x. y, z, T) = A( T, i) (9) 

is a function of T and i only. We store A( T, i) in relative location 167' + 
i of a block of 128 locations called the delta table. We can thus calculate 
the coordinates of the target atom by entering the delta table and using 
the addition 

h" = E+ A(T,i). (10) 

Similarly, the value of <f> a target atom is a function 0(7', *'), stored in 
relative location 167' + i of a block of 128 locations called the phi table; 
and this enables us to write <$>' into E' + 2 by a straight look-tip proce- 
dure. 
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The contents of the delta and phi tables are permanent settings com- 
puted from the lattice structure data cards before starting the first run. 
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